#PBAR_TestGaussExpFit.py
#
#
#  8/13/2013, John Kwong

def gauss_exp_function(x, a, x0, sigma, b, c ):
    return(a*np.exp(-(x-x0)**2/(2*sigma**2)) + b*exp(x*c))

def exp_function(x, b, c):
    return(a*exp(x*c))
    

# assumed load dataset from PBAR_ExamineBaseRecoverVer2.py
j = 0  # dataset index
k = 0 # detector number
i = 0 # time boundaries H: 0

y = dat[j][k][:,timeBoundariesList[i]:timeBoundariesList[i+1]].sum(1)
y = y / (timeBoundariesList[i+1] - timeBoundariesList[i])

# Set starting parameters
x = np.arange(256)
cut = (x > 25) & (x < 60)
startingParamGauss = [max(y[cut]), 40, 3]

# get starting parameters for the exponential
cut = (x > 60) & (x < 80)
pfit = np.polyfit(x[cut], log(y[cut]), 1)
startingParamExp = [exp(pfit[1]), pfit[0]]

startingParam = startingParamGauss + startingParamExp
cut = (x > 20) & (x < 70)
popt, pcov = curve_fit(gauss_exp_function, x[cut], y[cut], p0 = startingParam)

# 
plt.figure()
plt.grid()

plt.plot(x, y)
xarray = np.arange(80)
plt.plot(xarray, gauss_exp_function(xarray, popt[0], popt[1], popt[2], popt[3], popt[4] ))
plt.yscale('log')
